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ABSTRACT 

We have computed radiative equilibrium models for the gas in the circumstellar envelope sur- 
rounding the hot, classical Be star 7 Cassiopeia. This calculation is performed using a code that 
incorporates a number of im prov ements over previous treatments of the disk's thermal structure by 
iMillar fe Marlboroughl (| 19981) and| iJones. Sigut fc Marlboroug h ( 2004) ; most importantly, heating and 
cooling rates are computed with atomic models for H, He, CNO, Mg, Si, Ca, & Fe and their relevant 
ions. Thus, for the first time, the thermal structure of a Be disk is computed for a gas with a solar 
chemical composition as opposed to assuming a pure hydrogen envelope. We compare the predicted 
average disk temperature, the total energy loss in Ha, and the near-IR excess with observations and 
find that all can be accounted for by a disk that is in vertical hydrostatic equilibrium with a density 
in the equatorial plane of p(R) ~ 3 to 5 • 10~ n (i?/i?*) -2 ' 5 gem -3 . We also discuss the changes in the 
disk's thermal structure that result from the additional heating and cooling processes available to a 
gas with a solar chemical composition over those available to a pure hydrogen plasma. 

Subject headings: stars: circumstellar matter - stars: emission line, Be - stars: individual: 7 Cas 



1. INTRODUCTION 

Classical Be stars are non-supergiant B stars that pos- 
sess circumstellar material in the form of an equatorial 
disk. While the circumstellar disk is almost certainly a 
decretion disk of material from the star's atmosphere, 
the detailed mechanism that creates and maintains such 
a dis k remains unclear ([Porter fe Riv inus 200j| lOwockil 
l2004t ) . Rapid rotation of the central B star seems to play 
an important role, but there is still considerable debate a s 
to the extent ([Townsend et al.ll2004t iFremat et alJl2005[ ). 
Historically, the observational evidence for such circum- 
stellar material has been either spectroscopic or polari- 
metric in nature, and the accepted observational defini- 
tion of a Be star has been the appearance (at the current 
or previous epoch) of emission in the hydrogen B aimer 
lines. It has been recognized since the days of IStruvd 
(|1931| ) that recombination in a flattened disk could re- 
produce the range of spectroscopically observed Ha pro- 
files. In addition, the net (continuum) linear polariza- 
tion observed in Be stars is well explained by electron 
scat tering from non-spherically distributed circumstellar 
gas (|Covne fe Kruszewskil lT969: Wat ers fe Marlboroughl 

1199a) . 

Beginning with the resolution of </>Persei (B2 Vpe) 
at radio wavelengths with the Very Large Array by 
iDoughertv fe Taylor] (JT992) , interferometry has increas- 
ingly been used to spatially resolve circumstellar mate- 
rial. Be star disks have been resolved at radio, near-IR, 
and optical wavelengths, with these observations conclu- 
sively revealing the disk (lOuirrenbach et al.lll993l 119971 : 
iTvcner et al.l 120051 I2006D . T he observations are con- 
sistent with the suggestion of iPoeckert fe Marlboroughl 
(1978) that Be star disks are geometrically quite thin 
with opening angles of only a few degrees. 

Currently, optical interferometric observations require 
theoretical models of the emitting region in order to 



interpret the observed visibilities. Often observers fit 
simple models with free parameters to the data to de- 
scribe the disk emissivity. However, this simple pro- 
cedure can be considerably improved by using a de- 
tailed model for the thermal structure of the Be star 
disk. Such models naturally predict the emissivity 
and opacity of the gas r equired to produce theoreti- 
cal spectroscopic images (IMillar fe M arlborough [19981: 
Jones. Sigut fe Marlbor ough 2004; Carciofi & Bio rkmanl 



20061) . 



7 Cassiopeia (HD 5394; B0 IVe) is a interesting 
classical Be star which has a dens e, cool, equato- 
rial disk ([Millar fe Marlboroughl Il998t ). Thi s disk has 
been resolved with optic al interferometry (|Stee et alj 
119951: ITvcner et afl 120051 ). 7 Cas is likely the pri- 
mary in a binary syste m (Harman ec et al.l 120001 : 
iMiroshnichenko et al.ll2002D an d has unique X-ray char- 
acteristics ([Smith et al J 120041) . Although 7 Cas is of- 
ten quoted as "the prototypical" Be star, it has 
become clear that it possesses some unique char- 
acteristics. Nevertheless, it is a well-studied star 
making it an appropriate choice to test new codes 
and compare results with previously published work. 
In this current pa per, we extend the radia tive- 
equilibrium models of IMillar fe Marlboro ugh (1998]) and 
Uones. Sigut fe Marlboroughl (|2004f ) for the early Be star 
7 Cas to a disk models with a solar chemical composition. 

2. CALCULATIONS 

2.1. Overview 

The calculations in this paper were performed with a 
new code, bedisk, w hich is loosely based upon the cal- 
culational approach oflMillar fe Marlboroughl l[1998l ) and 
Uones. Sigut fe Marlboroughl (|2004l ). Models for the cir- 
cumstellar material were constructed assuming a den- 
sity distribution falling as an R~ n po wer-law in th e 
equatorial plane, following the models of Waters (1986), 
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ICote fe Watersi (|1987l ). and IWaters et all (|1987t ). By 
comparing to observations of the infrared excesses of 
Be stars, these authors found power-law density expo- 
nents in the range n ss 2.0 to 3.5. Perpendicular to the 
equatorial plane, it was assumed that the gas is in verti- 
cal, isothermal, hydrostatic equilibrium. Given the disk 
density distribution and the photoionizing radiation field 
from the central star, the equations of statistical equilib- 
rium were solved for the ionization state and level pop- 
ulations of H, He, CNO, Mg, Si, Ca, and Fe subject to 
the constraints of charge and particle conservation. Ra- 
diative transfer was handled via the escape probability 
approximation, and it was assumed that the dominant 
escape route for photons was perpendicular to the expo- 
nentially stratified disk. From this solution, the rates of 
energy gain and energy loss at each computational grid 
point were obtained, and the local temperatures were it- 
eratively adjusted to enforce radiative equilibrium. 

The main features and assumptions of this code are 
now discussed in detail. 

2.2. The bedisk Code 

The circumstellar disk was assumed to be axisymmet- 
ric about the star's rotation axis and symmetric about 
the mid-plane of the disk. In all that follows, R is used for 
the radial distance from the star's rotation axis and Z , 
for the perpendicular height above the equatorial plane. 
Thus the cylindrical co-ordinates of any disk location are 
(R,Z). If i?» denotes the stellar radius, then the calcu- 
lation domain is R* < R{ < i? m ax with i = 1 . . . n r , and 
< Zj < Z max (Ri) with j = 1 ... rig. Typically, n r is 
set to 60 with i? max /-R* = 50, and n z is set to 40. 

The code accepts a user-defined set of atomic mod- 
els which list the energy levels and the bound and free 
radiative and collisional transitions for each atom and 
ion to be included in the calculation. The set of atoms 
and ions and the total number of atomic levels and ra- 
diative transitions used for the current calculations are 
listed in Table [TJ For this initial work, the number of 
energy levels included for each atom and ion is similar t o 
the list of energy levels given bv lMoore fe Merr ill (1968). 
Sources for the required atomic data are given in Ap- 
pendix A. The abundances assumed for the various ele- 
ments can be found in Table [2] and are taken from the 
acc epted solar abundances of I Anders fe Grevessd (|1989l ) 
and lAnders fe Noelsl (|1993h . 

The number of atomic levels included for each 
atom/ion in Table [T] is fairly modest, although large 
enough to include most of the collisionally-excited lines 
seen in the optical and UV spectra of Be stars. While 
non-LTE solutions can be sensitive to t he number o f 
atomic levels included (see, for example, ISigutl (|1996l ) 
for a case-study of Cn in B stars), the computational 
time also increases sharply with the number of levels. 
Table [T] represents a compromise between realism and 
computational efficiency. Several techniques exist to 
group atomic levels into "super-l evels" (see, for example, 
iHubeny. Hummer. fe Lan3ll994[ ) which will allow future 
work to utilize more complete atomic models. 

The density structure of the disk is chosen in an ad- 
hoc manner. All calculations assume that the density 
drops as an R~ n power-law in the equatorial plane, and 
at each R, the gas is in vertical, isothermal hydrostatic 
equilibrium. To obtain the density above the equatorial 



TABLE 1 
Atomic Models 



z 


Atom 


Number of 








Levels 


Transitions 


\ 


H I 


15 




77 


\ 


H II 


1 







9 


He I 


1.3 




16 




He II 


1 







f; 


C I 


15 




2.3 


( : ; 


C ii 


10 




20 


5 


C in 


21 




45 


5 


C iv 


8 




15 


5 


C v 


1 







7 


N I 


23 




51 


7 


N ii 


15 




24 


7 


N in 


10 




20 


7 


N iv 


1 









O i 


15 




19 




O ii 


29 




80 




O in 


15 




21 


( ^ 


O iv 


10 




19 




O v 


1 







12 


A/To- t 


17 




28 


12 


Mg II 


12 




25 


12 


Mg in 


1 







11 


Si I 


37 




109 


14 


Si ii 


8 




1.3 


14 


Si in 


18 




27 


11 


Si iv 


11 




2.3 


14 


Siv 


1 







20 


Cai 


24 




42 


20 


Can 


11 




on 


20 


Cam 


1 







26 


Fel 


10 




90 


26 


Fen 


.39 




191 


26 


Fe III 


1 









TOTAL 


425 




1001 




TABLE 2 








Elemental Abundances 






Element 


Abundance 51 






H 




12.00 






He 




10.90 






C 




8.55 






N 




7.97 






O 

Mg 




8.87 
7.58 






Si 




7.55 






Ca 




6.36 






Fe 




7.51 





a The tabulated is abun- 
dance is log(Af/Af H ) + 12. 



plane at each radial distance, the user supplies a fixed 
set of density-drops from the equatorial plane, 

dj = \n{p{Z j )/p{Z 3 = 0)} (1) 

for j — l...nz- Here d\ — by definition, and 
d riz = —4. Then the Zj at which this density drop would 
occur, given the current value of Ri, is computed assum- 
ing vertical hydrostatic equilibrium with an isothermal 
temperature T , 

Z i= (( a/Rj \ 2 _\ h 

R, [{dj+a/Rj J ' W 
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Here a is given by 



GM, 



(3) 



where [i is the mean-molecular weight, w 0.5 for an ion- 
ized, pure hydrogen disk, and M is the mass of the central 
B-star. Thus the density at Zj) is given by 



(4) 



In this expression, p a , n, and T Q (through equation [3]) 
are the parameters which define the density structure of 
the disk. 

At each computational grid point, the photoionizing 
radiation field is required to evaluate the photoioniza- 
tion rates of all atoms/ions for the statistical equilibrium 
equations and to compute the photoionization heating 
rates for the radiative equilibrium solution. It is usual to 
divide this radiation field into a direct and diffuse con- 
tribution, 

J v = J° ir + J° if . (5) 

The direct contribution represents the radiation from the 
central star while the diffuse contribution arises from the 
disk. Note that despite this division, the only energy 
input into the circumstellar disk is assumed to be from 
the star itself so that ultimately, the energy in the diffuse 
field has its origin in direct radiation from the star. 

There is some evidence that 7 Cas has a binary 
companion with an o rb ital period of ~ 200 day s 
(jHarmanec et al.l 120001 : iMiroshnichenko et al.l l2002h . 
However, the exact orbit and the nature of the com- 
panion, including its spectral type and luminosity, are 
unknown. Rather than introduce further uncertain pa- 
rameters into the calculation (and invalidate the axisym- 
metric geometry), we shall assume that the energy input 
into the disk of 7 Cas from any potential companion is 
negligible. 

The direct component from the central star to the pho- 
toionizing radiation field at grid location (R4, Zj) is given 

by 



J v "(Ri, Zj 



(6) 



where G?f2 is an infinitesimal patch of solid angle centred 
around the direction h. The integral is over the visible 
stellar surface. Typically the surface is divided into a 
few hundred patches and the transfer equation is solved 
along a ray from the centre of each patch to the grid 
location 1 . The radiation field at the stel lar surface was 
taken from an LTE stellar atmosphere of iKuruci (|1993l ) 
which specified the mean intensity, J„(t„ = 0), at the top 
of the photosphere over a grid of 1221 frequencies. This 
was turned into the required intensity at each surface 
element by using the limb-darkening law 



V = !o {1 - au{\ - m) - Mi - m) 2 } 



(7) 



where the coeffic ients a u and b u wer e linea rly interpolated 
from Table V of I Wade fc Rucinskil (|1985l ). Here \i is the 

1 By choosing to solve the transfer equation along these rays, 
as opposed to simply applying exponential extinction of the stel- 
lar photospheric intensity, some contribution of the diffuse field is 
included. 



usual cosine of the surface viewing angle. Computing the 
mean intensity from this expression and setting it to the 
LTE model atmosphere prediction, we find that 



Io = 



■/„(t„ = 0) 
l-a,- (3/4) b v 



(8) 



While this procedure is approximate, and the exact 1^ 
for each fj, predicted by the LTE model could have been 
used, this approximation seems commensurate with oth- 
ers made in the construction of these models. We also 
note the use of a more physically realistic non-LTE, line- 
blanketed atmosphere for J v (t v = 0) would be a useful 
future improvement. 

In addition to limb darkening, gravity darkening in- 
duced by rapid stellar rotation can change the intensity 
distribution across the stellar disk and hence modify the 
direct component t o the photoionizing radiation field. 
iTvcner et al.l (|2005f ) interferometrically resolved 7 Cas's 
disk and estimated its inclination angle to the sky. They 
conclude the 7 Cas rotates at 0.7±0.1 of its critical veloc- 
ity. As 7 Cas does not seem to rotate particularly close to 
its critical velocity, we have not included gravity dark- 
ening (and the associated geometrical distortion of it's 
surface) into the calculation of the direct photoionizing 
radiation field. This point is further discussed in Sec- 
tion 13.11 where the adopted stellar parameters for 7 Cas 
are discussed. 

The simplest treatment for the diffuse field is to em- 
ploy the on-the-spot (OTS) approximation in which the 
recombination rate to level n of hydrogen is written as 



Rn,n( T n) — Rk.j 



(9) 



Here r ra is the optical depth at the continuum limit for 
photoionization from level n along a vertical ray to the 
nearest edge of the disk (this is consistent with our as- 
sumption that the dominant photon escape route is per- 
pendicular to the disk — see later discussion). Thus the 
principle assumption of the OTS approximation is that 
at high continuum optical depths, t„ 3> 1, recombina- 
tion to level n produces a photon that is locally absorbed 
within the same volume element, essentially undoing the 
recombination. Including the continuum optical depth 
dependence in equation © ensures that the OTS ap- 
proximation is used only when r„ > 1; the full recombi- 
nation coefficient is employed when the gas becomes op- 
tically thin in the continuum. We have applied to OTS 
approximation only to recombination to level n — 1 in 
hydrogen. The optical depths in the remaining continua 
are typically not large enough for a significant effect. We 
discuss a more complex, but still approximate, treatment 
for the diffuse field in section 13.51 

Given the photoionizing radiation field and current es- 
timates of the electron temperature and electron density 
at each grid location, we solve the statistical equilibrium 
equations to obtain the level populations for all atoms 
and ions. For atom k, and all of its associated ions, 
these are 

N k N k 

x»„ ( 10 ) 



where is the number of atomic levels 



for i = 1 

included for the k th atom and its ionization stages. Note 
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that these equations must be supplemented by a particle 
conservation equation of the form 



(11) 



for each atomic species. The elemental abundance, Ak, 
can be found from Table [U 

In the case of bound-bound transitions i — ► j, i < j, 
the rates have the simple form 



Rij = n e qij(T e ) (12) 



and 



Rji = Aji P csc (t z ) + n c qji(T e ) . (13) 

Here the factors qij(T e ) and qji(T e ) represent collisional 
excitation and de-excitation respectively and are propor- 
tional to the Maxwellian-averaged collision strength for 
i — ► j transition. Aji is the usual Einstein transition 
probability for spontaneous emission. This form of the 
statistical equilibrium equations handles radiative trans- 
fer in the line via the escape-probability approximation. 
The escape probability for each grid location was ob- 
tained by computing the line-centre optical depth to the 
nearest vertical edge of the disk (denoted r z ) and then us- 
ing this optical depth to estimate the static, single-flight 
escape probability assuming complete redistribution in 
the spectral line. 

We have assumed that the dominant loss route for pho- 
tons is perpendicular to the disk because of the exponen- 
tial density stratification implied by vertical hydrostatic 
equilibrium. As it is reasonable to assume that the main 
motion of the disk gas is Keplerian rotation about the 
central star, the Doppler shifts experienced by photons 
escaping roughly perpendicular to the disk will be small, 
and it is appropriate to approximate the escape proba- 
bility, -Pesc(Tz), by a static, single-flight escape probabil- 
ity, as opposed to employing the Sobolev approximation. 
This is consistent with the definition of t z as the optical 
along a ray perpendicular to the disk to the nearest edge; 
t z is not defined in terms of a local velocity gradient as 
in the Sobolev approximation. 

The form of the static, single-flight escape probability 
appropriate for complete-redistribution over a Doppler 
profile is 

Pcsc(t) = ; , 1 — , (14) 



4r(ln(r/^F)) 1/2 ' 



where r > 1 (|Mihalasl 119781 [C anno 5 H985h . As this 
result is an asymptotic expansion (essentially in the limit 
that the scale of variation of the line source function is 
large compare d to the width of the scattering kernel - see 
ICannonlll98 5). an ad- hoc correction is needed to handle 
r <C 1 . Here we have adopted the suggestion of iTielensI 
(|2005ft where the escape probability is taken to be 



P csc (t) 



-2.34t 



4.68t 



(15) 



for r < 7. This function is continuous with the asymp- 
totic result at t = 7 and tends to the limit of 1/2 as 
r — > 0, a reasonable result if the collisional destruction 
probability of the line photon upon scattering is not too 
small. 



In the case i — > j = k is a bound-free transition, the 
photoionization and recombination (spontaneous plus 
stimulated) rates are given by 



dv 
hv 



and 



= 4?r 



Ju 



dp 



(16) 



(17) 



Here (rii/n K )* is the LTE population ratio found from 
the Saha-Boltzmann eq uation and it i s proportional to 
the electron density, n e (|Mihalas l ll978h . Dielectronic re- 
combination and autoionization are included by retaining 
the full resonance structure of the photoionization cross 
section ai K {v). 

Given the solution for the atomic level populations, a 
new estimate for the electron density can be made by en- 
forcing charge conservation, and the rates of heating and 
cooling for the various atomic processes can then be com- 
puted. Heating includes photoionization and collisional 
de-excitation while cooling includes radiative recombina- 
tion and collisional excitation. Deta iled express i ons fo r 
all of these processes can be found inlpsterbrockl (j 1989f> . 
However, in contrast to lOsterhrockl (|1989l ). transitions 
(and the implied cooling) due to radiative recombina- 
tion were computed explicitly via equation (|17[) for each 
atomic level; total recombination co-efficients (summed 
over n) were not used. 

Heating due to viscous dissipation in a K eplerian disk 
was also included (|Lee. Saio. fe Osaklll991[ ) but was al- 
ways found to be negligible. 

Net cooling due to free-free emission (in the fields 
of Hi and Hen) w a s incl uded via the expression of 
Rvbicki fc Lightmanl (|1979h . modified as suggested by 
Netzerl ( 1990T ) to account for the reduction in the free- 
free cooling rate as the gas becomes optically thick to 
free-free radiation, 



L ff = 1.4 • 10" 



T K n, 



E 

i=H,He 



9ff< 



K /kT s 



(18) 

He re Zj = 1 and the frequency cut-off, i/ max , suggested 
by iNetzerl (|1990t ). is the smallest frequency for which 
the optical depth to the nearest edge of the disk ex- 
ceeds one. The Gaunt factor, <?//, which varies slowly 
with temperature, was set to a constant value of 1.2 
which iRvbicki fc Lightmanl ()1979f ) indicate will approxi- 
mate the exact result to within 20%. 

To find the equilibrium kinetic temperature, T e , at 
each grid location, heating and cooling were balanced 
by searching for a zero in the net-cooling rate, r]c(T e ). 
The root was initially located via bisection and then re- 
fined with the secant method. In the rare case of multiple 
roots, the stable one satisfying drjc/dT e > was chosen. 

The overall flow of the calculation is to start at the 
inner boundary of the disk, closest to the star at i = 1. 
Solutions proceed downward in Z, from the top of the 
disk (j = n z ) to the equatorial plane (j — 1). This allows 
the optical depths back to the star and to the nearest 
edge of the disk to be kept current with the solution 
level populations. 
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3. COMPUTATIONS 

In this section, we first compare the predictions of our 
code with known results for 7 Cas. Next we explore a 
wide range of disk parameters for 7 Cas and investigate 
the effect on the temperature structure of the disk of us- 
ing a solar chemical composition for the gas. We then 
examine the energy loss in the spectral lines included in 
the models and compute the near infrared spectral en- 
ergy distribution. Finally, we examine the computation 
of the diffuse photoionizing radiation field generated by 
the disk itself in order to evaluate the use of the OTS 
approximation for most of the models computed in this 
work. 

3.1. Comparison with Millar & Marlborough 

iMillar &; Marlborough! (l998) constructed a pure- 
hydrogen radiative equilibr ium model for 7 Cas. 
IMillar fc Marlborough] (|l999h (MM, hereafter) extended 
this work to include the OTS approximation for the dif- 
fuse radiation field, and it is this temperature distribu- 
tion which we have chosen for comparison. We adopt 
the same stellar parameters for 7 Cas as MM, which are 
reproduced in Table [3] As we use these fundamental pa- 
rameters for all of the calculations in this work, some 
additional comment is in order, particular ly concerning 
the ad opted stellar effective temperature. iFremat et al.1 
( 2005) investigate of the effect of rapid rotation on funda- 
mental parameter determinations of B stars. They find, 
by fitting the line spectrum between 4250 and 4500 A, 
"apparent stellar parameters" for 7 Cas (those obtained 
by a best fit classical, plane-parallel model atmosphere) 
of T e ff = 26, 400 K and log(g) = 3.8 which are close to the 
parameters adopted in this work. Nevertheless, they do 
find significant effects of rotation in their best-fit rotat- 
ing models which have a parent non-rotating counterpart 
T c s (see their paper for details) of sa 30, 000 K. This 
result suggests that acc ounting for th e rotat ion of 7 Cas 
m a manner following IFremat et ail (|2005l ) (but using 
non-LTE stellar atmospheres) would be a useful future 
improvement in the computation of the direct stellar con- 
tribution to the photoionizing radiation field. 

The fixed density struct ure for the 7 C a s disk adopted 
by MM is described by iMarlboroughl (|1969f ) and is 
slightly different from the model described in section l2~2l 
MM assume that the disk is in isothermal, hydrostatic 
equilibrium at only one radial distance from the central 
star and that the radial drop-off in the equatorial density 
follows from an assumed (radial) outflow velocity law and 
the equation of continuity. We have simply adopted the 
(Ri, Zj) grid of MM (which is 24 by 20) and their total 
density at each grid point as input to BEDISK. We have 
used a 5-level hydrogen atom plus continuum for this 
comparison and have also used the OTS approximation 
(as described previously) for the diffuse field. Despite 
this, there are still some significant differences between 
the calculations: our approach uses the optical depths 
in the OTS approximation and in the line escape prob- 
abilities as opposed to the various cases of MM. We use 
a newer ATLAS stellar atmosphere to predict the pho- 
toionizing radiation field from the star and use many 
more rays from each grid point back to the star. No 
attempt was made to use identical atomic data: MM 
included collisional transitions for only transitions n to 



TABLE 3 
Stellar Parameters for 7 Cas. 



Parameter 


Value 


Unit 


Spectral Class 


B0 IVe 




Radius 


10.0 




Mass 


17.0 


Mq 


Luminosity 


3.4 10 4 


Lq 


T G ff 


25,000 


K 


l°g(s) 


3.50 


cms" 


Distance" 


1881^8 


pc 



a Dista nce from the Hipparco s cat- 
alogue (Pcrryman et al. 1 99fl) . 

n ± 1 in hydrogen whereas we have included all colli- 
sional rates. Nevertheless, despite these differences, the 
comparison is a useful check. 

bedisk predicts a density-weighted average tempera- 
ture, defined as 

%=^f [ TpdV, (19) 

M J Disk 

(where M is the total mass of the disk), of 11300 K. 
This is to be compared to the 14 500 K quoted by MM 
(in their Table 1). This significant difference is simply 
one of definition: the density-weighted average quoted by 
MM is actually defined as (^ pij Tij ) / Pij where the 
sum is over all of the grid points in the calculation. Com- 
puting this quantity for the current bedisk model yields 
13 900 K which agrees with the MM result to within 5%. 

Figure[T]compares the ratio of the bedisk temperature 
to the MM temperature throughout the entire circum- 
stellar disk. Agreement is generally good; bedisk tends 
to be somewhat cooler near the equatorial plane in the 
inner portion of the disk, while somewhat hotter towards 
the upper edge. However, 80% of the grid points agree 
to within ±20%. The largest differences tend to occur 
along with upper edge of the envelope where the opti- 
cal depths (to the nearest vertical edge of the disk) are 
most rapidly changing. The treatment of the escape of 
line radiation, as noted above, and particularly of how 
the OTS approximation was implemented (MM applied 
OTS to the whole disk as opposed to including an optical 
depth dependence as in equation [5]) , are likely the origin 
of the more significant differences. 

3.2. Effect of Adding Metals on the Thermal Structure 

Table 2] gives the disk parameters for 16 disk models 
computed to compare with observations of 7 Cas. These 
models span a range of nearly two orders of magnitude 
in density (as obtained by varying p Q from 2.5 ■ 10~ 12 to 
1.0- 10~ 10 gmcm~ 3 ) with two values assumed for the ra- 
dial drop-off of the density in the equatorial plane, R~ 2 5 
and R- 3 - 5 . All models assumed T Q = 13 500 K for the 
isothermal temperature which sets the vertical density 
scale-height via Eq. [3J Also given in the table are the 
predicted density- weighted temperatures, the disk emis- 
sion measures (in cm -3 ), defined as 

EM = f n%dV, (20) 

./Disk 

and the predicted total Ha luminosities (in ergss -1 ). 
Section 13.31 discusses how the Ha luminosity was com- 
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Fig. 1. — Ratio of the BEDISK temperature to the MM tempera- 
ture for 7 Cas. 



TABLE 4 
Model Disk Parameters. 



n 


Po 






log EM 


log Ha 




g cm" 


3 


K 


cm ° 


erg s 


2.5 


2.510" 


12 


14060 


59.25 


33.43 


2.5 


5.010" 


12 


12870 


59.84 


33.79 


2.5 


7.510" 


12 


12420 


60.16 


33.99 


2.5 


1.010" 


11 


12170 


60.40 


34.12 


2.5 


2.510" 


11 


11040 


61.10 


34.48 


2.5 


5.010" 


11 


9420 


61.50 


34.69 


2.5 


7.510" 


11 


8590 


61.57 


34.77 


2.5 


1.010" 


10 


8140 


61.62 


34.80 


3.5 


2.510" 


12 


13990 


58.88 


32.90 


3.5 


5.010" 


12 


13740 


59.48 


33.17 


.3.5 


7.510" 


12 


13500 


59.81 


33.31 


3.5 


1.010" 


11 


13290 


60.05 


33.40 


3.5 


2.510" 


11 


12080 


60.82 


33.68 


3.5 


5.010" 


11 


11050 


61.35 


33.85 


3.5 


7.510" 


11 


10560 


61.48 


33.95 


3.5 


1.010" 


10 


10240 


61.55 


34.01 



puted. The parameters of the central star were again 
those of Table [3 

The density-weighted temperatures predicted by the 
models are plotted in Figure [2] as a function of p . These 
results are compared to the predicted density-weighted 
temperatures for a set of pure hydrogen disks with identi- 
cal physical parameters. Also shown in the Figure is the 
observed disk temperature for 7 Cas of 9500 ± 1000 K as 
found bv lHonv et al.l (|200 Q:) by fitting the IR Humphrey's 
bound-free jump at 3.4 /im. As expected, denser disks 
predict lower density- weighted temperatures. 

The R~ 2 5 models are consistent with the iHonv et al.l 
(2000) result for densities in the range of 3 to 8 • 
10~ 11 gc m~ 3 . This agrees well with the density esti- 
mated bv lHonv et all (120001) using their observations and 
the disk models of lWaterd (119861). The i ? ~ 35 m odels are 
only just consistent with the lHonv et all (|2000D result for 
largest densities considered, p sa 10~ n gcm~ 3 

The predicted temperature trend as a function of p a 
has an interesting dependence on the metallicity of the 
gas. For low density disks, the solar composition disks 
are considerably cooler than the pure hydrogen models 
by 1-2000 K. However, the difference decreases for higher 
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Fig. 2. — Density- weighted disk temperatures for the models of 
Table H 

disk densities, p > 2- 10~ n gem 3 . Indeed, for the i?~ 2,5 
models, the higher density solar and pure hydrogen disks 
predict nearly the same density- weighted temperatures. 
This behaviour can be understood in terms of the heat- 
ing and cooling avenues introduced by metals. Metals 
can act to cool the gas due to the escape of collisionally- 
excited line radiation. However, metals can also help to 
heat the gas via photoionization. If the optical depths in 
the hydrogen continua (excluding the Lyman continuum) 
are low, then the additional heating provided by the pho- 
toionization of metals is negligible in comparison to hy- 
drogen. In this case, it is the cooling due to collisionally- 
excited line radiation that dominates. In a low-density 
disk, line cooling is further enhanced by the small optical 
depths (and hence high escape-probabilities) in the lines. 

At high densities, however, the optical depths in the 
hydrogen bound-free continua are much larger; pho- 
toionization heating due to metals can then become 
important, particularly as many abundant metals have 
bound-free thresholds in the short-wavelength region of 
the Balmer continuum which is near the photospheric 
flux maximum in B stars. In this case, metals add 
both heating and cooling and the net result is a very 
similar density-weighted temperature to the case of a 
pure hydrogen plasma. The se trends help explain why 
(Millar & Marlborough 1998]) where able to obtain a rea- 
sonable density-weighted temperature for 7 Cas despite 
using a pure hydrogen envelope. 

Two-dimensional temperature distributions in the disk 
for several different density i?~ 2 5 models are shown in 
Figured] The inner portion of the disk for R/R* < 5 is 
expanded in each case for clarity. These figures clearly 
show the development of a cool region near and in the 
equatorial plane for the higher density disks. These 
denser disks have fairly strong vertical temperature gra- 
dients perpendicular to the equatorial plane. Given this, 
it would seem prudent to re-integrate the equation of hy- 
drostatic equilibrium at each radial distance, accounting 
for the vertical variation of the gas temperature, and then 
to iterate the pressure structure along with the thermal 
solution to produce a disk that is in both radiative and 
(vertical) hydrostatic equil ibrium; we shall present suc h 
models in a future work (|McGill. Sigut fc Jonesl 12006) . 
However, as the focus of the present work is on the ther- 
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Fig. 3. — Circumstellar disk temperature as predicted by three 
ij -2,5 models. The top panel has p a = T 10 — 12 gem -3 , the middle 
panel, p a = 5 • 10 — 12 , and the bottom panel, p = 5 ■ 10 -11 . Each 
panel on the right enlarges the inner portion of the disk from the 
panel on its left. 
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Fig. 4. — Temperature ratio between a full solar composition 
disk to that of a pure hydrogen disk. Shown in the bottom panel 
is a model with p = 5 ■ 10 — 12 gem -3 , and in the top panel, one 
with p = 5 ■ 10 -11 gcm -3 . Both models had an R~ density 
drop-off in the equatorial plane. Each panel on the right enlarges 
the inner portion of the disk from the panel on its left. 



mal structure of the disk, it is convenient to have a fixed 
density structure so that the thermal effect of the gas 
metallicity can be unambiguously seen. 

Figure [4] presents a detailed comparison of the tem- 
perature structure predicted by a gas with a realistic so- 
lar composition to that of a pure hydrogen model. The 
figure compares the temperature ratio T Solar /T PuTeii f or 
two densities, p = 5- 10 -11 and p Q — 5- 10~ gem -3 . In 
both cases, the optically thin gas far above the equatorial 
plane is cooler with the inclusion of metals. However, in 
the lower density model, p = 5- 10 -12 gem -3 , the gas in 
the equatorial plane near the star is hotter for the solar 
composition gas out to R/R* ~ 6. In the higher density 
model, the gas near the equatorial plane has nearly the 
same temperature in the solar and pure hydrogen mod- 
els. These trends are consistent with the effects noted in 
the discussion of the density-weighted average disk tem- 
peratures. 




■1U -10.6 

logP D (g/cm 3 ] 

Fig. 5. — The predicted energy loss in Ha from Eq. 1211 (in 
ergss - 1 ) for the models of Table|4] Shown as the shaded rectangle 



Mazzlil tl989ft and 



is the observed r ange of values from IKastner - 
llStee et al.H1995l) . 



3.3. Line Luminosities 

Figure [5] plots the total energy lost in Ha (in ergss -1 ) 
by the 16 models of Table 01 The line luminosity was 
found by integrating the flux divergence, in the escape 
probability approximation 2 , over the volume of the disk, 



Disk 



7>j Pesc ( 



i)dV. 



(21) 



For Ha, i is level n = 2 of H I and j is level n = 3. Also 
shown in the F igure is the Ha lum i nosity obse rvationally 
determ ined by iKastner &; Mazzlil (|1989f ) and IStee et al.1 
(1995). The Ha luminosity from 7Cas is known to be 
variable, but the cited values are typical of the current 
epoch. There is good agreement with the observed lumi- 
nosity for R~ 2 b models with p a between 2.5 • 1CP 11 and 
10 -11 gcm -3 . This result is consistent with the models 
that b est fit the observed disk temperature of lHonv et al.l 
(2000). However, the i? -35 models seem inconsistent 
with the total energy loss in Ha for the range of disk 
densities considered. 

Figure [5] shows energy escaping in all of the included 
radiative transitions for the R~ 2 5 model with p Q = 
5 • 10~ n gem -3 . The fluxes are again found by integrat- 
ing Eq. [21] over the disk. It is important to keep in mind 
that Figure[6]is not a spectrum (which would be obtained 
by integrating the transfer equation along a series of rays 
through the computational domain); it is simply a plot 
of the energy loss per second in each line acting to cool 
the gas. Nevertheless, it gives a good indication of the 
expected strong emission lines in the disk spectrum. For 
this particular model, 94% of the energy loss is provided 
by the lines of Hi. Contributing at the level of « 1% 
percent are Fen, C11, Mgn, and Hei. Next to the lines 
of Hi, the largest energy losses are in the resonance lines 
A 1333.6 A line of C 11 (2s 2 2p 2 P° - 2s 2p 2 2 D) and the h 
and k lines of Mg 11 near A 2800 A. Although it does not 
possess a single strong line, cooling due to Fen domi- 
nates over all of the metals due to its rich spectrum with 
many collisionally-excited lines. It should be noted that 
Figure [5] and the percentage contributons cited above 

2 In the escape probability approximation, the net radiative 
bracket is replaced by the (single-flight) escape probability. 
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The energy loss in each radiative transition predicted by 



Fig. 

the R~ 2 5 disk model with 



Po 



:5 - 10" 



Only transitions 



with an energy loss above 10 30 ergs s — 1 are shown. The symbols 
representing the various elements are as indicated in the legend. 
The ionization stage is coded by colour: neutral (black), first ion 
(blue), second ion (red), third ion (green). 

represent a global picture of energy loss integrated over 
the entire disk. The impact of the heating and cooling 
contributions of metals can be much larger at individual 
grid locations as demonstrated by Figure |U 

3.4. Predicted IR spectral energy distributions 

In this section, we consider the predicted IR contin- 
uum energy distribution as emitted perpendicular to the 
disk 3 . Such models have a zero inclination angle between 
the rotation axis and the observer's line of sight. This 
spectrum is easily found by solving the radiative trans- 
fer equation vertically through the disk at each If 
Ii = I + i v (Ri) is the emergent intensity for the annulus 
of area A{ = ""(iZf+x/a — then the SED seen by 

the observer will be 



- Star+Disk 



11 n 

E 

i=l 



(22) 



Here i?* is the radius of the star, /* is the specific in- 
tensity corresponding to the stellar surface, and Ii is the 
specific intensity of the i th disk annulus. It is well known 
that at IR wavelengths, Be stars possess an excess of radi- 
ation over that predicted by an appr opriate stellar photo- 
sphere model (jCote fc Waterslll987t ). This excess comes 
mainly from free- free emission in the ionized disk 4 , with 
a small contribution from free-bound emission shortward 
(in wavelength) of the ionization edges. 

The near-IR spectral energy distributions are shown in 
Figure for the R- 2 b models of Table |U Plotted is the 
(monochromatic) IR excess expressed in magnitudes, 



2.5 log 



j^Star-\-Disk 
J^Star 



(23) 



The discontinuous jumps at wavelengths less than 5 pxa 
in this figure, particularly for higher densities, repre- 
sent the hydrogen free-bound continua for recombina- 
tion to n = 4 (1.4/im) and n — 5 (2.2/im). Unlike 

3 We shall present detailed synthesis for the infrared hydrogen 
disk spectrum for arbitrary inclination angles in a future work. 

4 Note that the exact expression for the free-free Gaunt factor 
was used for this calculation. 



(plane-parallel) stellar photospheres which have bound- 
free edges in absorption (reflecting the inward increase in 
temperature), circumstellar material exhibits the bound- 
free edges in emission, reflecting the increase in the gas 
emissivity due to recombination. Thus the figure indi- 
cates a small contribution from free-bound emission to 
the infrared excess at these wavelengths. We also note 
that as the disk density is increased, the infrared excess 
steepens most strongly between the wavelengths of 1 and 
5 jum. 

To compare with observations, we first show the 
12 and 25 11m I RAS IR excesses for 7 Cas found by 
ICote fc Waters! (|1987l ). The model that best matches 
the IRAS IR excess at these wavelengths has a slightly 
smaller density, p Q ss 10 -11 gem -3 , than the model that 
best matches the observed average disk temperature and 
energy loss in Ha, p Q ss 3- 10 -11 gem" 3 (Figure^]). How- 
ever, there is good ev idence that 7 Cas is seen at an in- 
clination of i >« 55° (|Tvcner et al.1 12006) . as opposed to 
being viewed pole-on (i = 0°) as assumed by the models. 
Hence the models have an effective emitting area that is 
too large and the IR excess is likely overestimated. In 
addition, the spectra for non-zero inclination angles will 
reflect the contribution of a different set of rays passing 
through the disk; the different physical conditions along 
these rays will lead to differences in the gas opacity and 
emissivity and hence a different predicted intensity along 
each ray. 

We also show in Figure [7] the 2.5 to 11.6 /im Infrared 
Space Observatory (ISO) spectrophotometry 5 for 7 Cas. 
In order to compare the ISO fluxes at the Earth to our 
photospheric model (to set the reference flux to extract 
the excess), we requ ire the angular diameter of the star. 
iTvcner et al.1 (|2005h cite the major axis of the Ha emit- 
ting region of 7 Cas of 3.67 ± 0.09 milli-arcseconds, but 
this has a large contribution from the circumstellar disk. 
To resolve this problem, we have chosen an angular diam- 
eter such that the ISO data reproduces the 12 /im IRAS 
excess. Combining this stellar angular diameter for 7 Cas 
with the range of Hipparcos distances listed in Table [3l 
we find a required radius for 7 Cas of between 6.8 and 
8.4 solar radii. While this result in disagreement with 
the 10 solar radii listed in Table 02 the larger value is 
within the plausible error range for 7 Cas's radius. With 
this normalization, the ISO observations match reason- 
ably well with the p ~ 10" n gcm~ 3 model prediction 
over the considered wavelength range. The fit is actu- 
ally worse closer to the normalization point; the shorter 
wavelength excess, between 2.5 < A < 5 ^m, is quite well 
reproduced. Similar caveats to those ending the previ- 
ous paragraph apply here: a detailed comparison will be 
made with a model that correctly accounts for the view- 
ing inclination of 7 Cas in a future work. 

3.5. An Approximate Treatment of the Diffuse 
Radiation Field 

While the OTS approximation is simplest treatment 
of the diffuse radiation field, Figure [3] shows that in the 
densest disks (p « 5 ■ 10 _11 gcm -3 ), grid locations near 
the equatorial plane can become quite cool due to large 

5 We have extracted this data from the ISO on-line archive. The 
spectrophotometry is from the PHT-40 instrument for observing 
series TDT 76803401. 
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FlG. 7. — The IR excess expressed in magnitudes predicted by 
the R- 2 - 5 mod els of Table [j] The sq uare data points at 12 and 
25 fim are from [Cote & Waters (1987). The open circles are the 
ISO spectrophotometry. 




R/R„ 

FlG. 8. — The effect of adding the diffuse contribution of Eg. 1251 
to the temperatures predicted for the R~ 2 5 disk model with p = 
5 ■ 10 -11 gc m~ 3 . Plotted is the temperature ratio of a calculation 
including Eg. 1251 to one omitting this contribution. 



optical depths back to the central star. However, exam- 
ining the thermal structure in the Z-direction at such a 
location typically shows that at heights above and be- 
low the equatorial plane, the gas can still be quite hot 
as it can be directly illuminated by at least a portion of 
the star. One might expect that the radiation emitted 
from these hot "sheaths" might be an important source 
of s econdary photoionizing rad iation for the equatorial 
gas ()Carciofi fc Biorkma n 2006). 

In principle, the calculation of J^ lf in Eq. [5] requires a 
solution of the transfer equation in the 2D cylindrical ge- 
ometry 6 . However, to estimate the potential effect of this 
diffuse field on the current work, we have tried a simpler 
and approximate approach. We solve, at each Ri, the 
radiative transfer equation perpendicular to the disk in 
the Z-direction using th e method of short-characteristics 
(01 son fc Kunaszlfl98l . The mean intensity along this 
ray is then 

MZj) = \ {I^+iAZj) + I,=-iAZ 3 )} ■ (24) 

For this solution, zero incident radiation is assumed per- 
pendicular to the disk, We then estimate the diffuse con- 
tribution to the photoionizing radiation field at each Zj 
by assigning 

J™{Z 3 ) = Wdif Ju{Zj), (25) 

Here W<m represents an ad-hoc dilution factor, between 
zero and one, for the perpendicular rays. The expecta- 
tion is that J^* lf will potentially be most important in 
the cool equatorial regions that develop for higher den- 
sity disks such as the one illustrated in the bottom panel 
of Figure E3 In such a cool obscured region, the hot 
sheaths above and below cover approximately 1/2 of the 
available 47r steradians so a dilution factor of Wdif = 0.5 
is a reasonable approximation. In regions of the disk 
where the central star is not obscured, the choice of Wdif 
is irrelevant as J Dlr 3> J Dlf . For this reason, we have 
taken Wdif to be 1/2 at all grid locations in the disk. 

6 A direct implementati on of this approach can be found in 
Carciofi & Biorkman (2006) who use a Monte Carlo approach to 
radiative transfer to estimate the diffuse radiation field in the en- 
velope of a somewhat later-type Be star. 



Figure [H] shows the effect of adding such an approx- 
imate treatment of J to the R~ 2 5 disk model with 
po = 5 • 10~ n gcm~ 3 , the model which best matched 
the observed disk temperature and the total Ha lumi- 
nosity. With the diffuse field added, the amount of ad- 
ditional heating in the equatorial plane is fairly small, 
with only a 10 — 20% increase in the temperature there. 
Hence, we feel that the previously computed models are 
reliable despite the approximate treatment of the diffuse 
field through the OTS approximation. 

Note that this figure is not meant to imply that the 
effect of the diffuse field is negligible; in both runs, the on- 
the-stop (OTS) approximation, which approximates the 
local trapping of radiation with A < 912 A and hence the 
diffuse photoionizing radiation field for A < 912 A was 
employed following Eq. [9] If the OTS approximation is 
turned off, there i s considerable cooling of the e quatorial 
regions (see also iMillar fc Marlborough! 1 19991) and this 
would mask the effect (illumination of these same regions 
from above and below) that we were trying to investigate. 

4. CONCLUSIONS 

For the first time, we have constructed radiative equi- 
librium models for the disk of 7 Cas using a gas with 
a solar chemical composition. We find that a model 
with a power-law equatorial gas density of p(R) ~ 3 to 
5 • 10 _11 (i?/i?*) _2 ' 5 gem -3 , which is in vertical (isother- 
mal) hydrostatic equilibrium, can reproduce several 
overall observed disk properties, including its (density- 
weighted) average temperature, the energy loss in Ha, 
and the near infrared flux excess. 

Also in this work, we have investigated the differ- 
ences between our solar composition radiative equilib- 
rium models and the corresponding set of pure hydrogen 
models with the same disk parameters. The effect of 
the additional heating and cooling provided by elements 
heavier than hydrogen on global measures of the disk 
structure, such as the density-weighted average disk tem- 
perature or the total energy loss in Ha, depend strongly 

7 A s ubtle point is whether a model including J Dlf according to 
Eq. 1251 should also employ the OTS approximation. In principle, 
employing both over counts the local radiation field in each volume 
element. 
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on the parameter p Q (the assumed equatorial density near 
the stellar surface) which sets the overall density of the 
disk. For low disk densities, p < 10~ n gcm~ 3 , the 
density-weighted averaged disk temperature is generally 
1-2000 K cooler with a solar composition, due mostly 
to enhanced collisionally-excited line cooling. At higher 
densities, absorption of ionizing radiation in the bound- 
free continua of elements heavier than hydrogen provides 
additional heating which offsets somewhat the additional 
line cooling. 

In examining the detailed temperature at each position 
withing the disk, differences of up to w ± 40% are found 
in comparing the solar composition models to the pure 
hydrogen models. While the global energy loss in Ha in 
the solar and pure hydrogen models is similar, these sig- 
nificant differences in the disk temperature distribution 
may manifest themselves in the detailed spectra. The 
next step in the comparison between these two sets of 
models is to present predicted spectra over a wide range 
of wavelengths for a wide range of viewing inclinations. 

The bedisk code represents a compromise between 
computational efficiency and realism. The code's most 
notable approximations are the treatment of the disk dif- 
fuse radiation field, the use of (first-order) static, escape 
probabilities for the line radiative transfer, the assump- 
tion of a spherical, non-rotating, central star, and the 
somewhat limited atomic models. Nevertheless, these ap- 
proximations allow the code to efficiently explore a wide 



region of parameter space for the thermal structure of the 
circumstcllar disks surrounding hot stars. Future work 
will proceed along two fronts: we will use the current 
version of bedisk to compute a series of Be disk thermal 
models that will serve as the basis for detailed spectral 
synthesis to produce line profiles, interferometric visibil- 
ities, and continuum polarization signatures which can 
be directly compared to observations. Along the second 
front, we shall improve some of the basic assumptions of 
the code. Most notably, we will enforce a vertical hydro- 
static equilibrium density structure consistent with the 
thermal solution, allow for potential gravitational dark- 
ening (and geometric distortion) of the central star by 
rapid rotation, and account for the diffuse photoionizing 
radiation field from the disk by a direct solution of the 
2-D radiative transfer problem. 
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APPENDIX 
THE ATOMIC DATA 

Energy levels for all atoms and ions were adopted from the NIST Atomic Spectra Database 8 . Radiative transition 
probabilities and photoionization cross sections were adopted from the Opacity Project Database TOPbase 9 . The 
photoionization cross sections, which include complex resonance features due to autoionization, were smoothed by 
convolution with a Gaussian down to the resolution of the ATLAS frequ ency grid . The r mally-averaged collisio n 
stre ngths for the electron impact ex citation of hydrogen were adopted from ICallawavi (|1994h . lAggarwal et al.l (|1991h . 
and lChang. Avrett &: Loeserl (|1991|). Thermally-averaged collisional strengths for the remaining atoms and ions were 
adopted from the compilation o f lPradhan fcPenel (fl995h where available. The Gaunt factor approximation was used 
for the remaining dipole-allo wed transitions. The at omic data for iron and its ions were adopted from the extensive 
model atoms constructed bv lSigut fc Pradhanl |2003) and lSigut. Pradhan k. Nahari (|2004h . 
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